LTER Pulses WaveletComp R Package Tutorial
KM Hall
2022-09-14
Refer to guided tour document
Basic background for understanding a wavelet transform of a time series
library(WaveletComp)
library(dplyr)
library(matrixStats)
library(tidyr)
library(ggplot2)
library(lubridate)
library(readr)Load the data
seg <- read_csv("~/Documents/LTER_Pulse_Time_Series_WG/SEV_daily_flux/US-Seg_daily_aflx.csv") %>%
mutate(date = mdy(Date)) %>% # WaveletComp expects date var named "date"
select(-Date) %>%
select(date, everything()) %>%
arrange(date) # make sure data are properly sorted by date## Rows: 5114 Columns: 23
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Date
## dbl (22): GPP_int, GPP_dayint, RECO_int, RECO_dayint, RECO_nightint, NEE_int...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.str(seg)## tibble [5,114 × 23] (S3: tbl_df/tbl/data.frame)
## $ date : Date[1:5114], format: "2007-01-01" "2007-01-02" ...
## $ GPP_int : num [1:5114] 0 0.00395 0.0362 0.25915 0.46783 ...
## $ GPP_dayint : num [1:5114] 0 0.00228 0.0362 0.13995 0.24107 ...
## $ RECO_int : num [1:5114] 0.495 0.61 0.7 0.592 0.514 ...
## $ RECO_dayint : num [1:5114] 0.247 0.262 0.234 0.273 0.268 ...
## $ RECO_nightint: num [1:5114] 0.249 0.348 0.465 0.32 0.246 ...
## $ NEE_int : num [1:5114] 0.495 0.606 0.664 0.333 0.046 ...
## $ NEE_dayint : num [1:5114] 0.2467 0.2594 0.1982 0.1328 0.0272 ...
## $ NEE_nightint : num [1:5114] 0.2486 0.3464 0.4654 0.2005 0.0188 ...
## $ ET_int : num [1:5114] 0.368 0.25 0.255 0.458 0.864 ...
## $ ET_dayint : num [1:5114] 0.359 0.262 0.305 0.425 0.794 ...
## $ P_int : num [1:5114] 0.614 1.024 0 0 1.024 ...
## $ TA_avg : num [1:5114] -8.53 -5.98 -7.56 -1.98 3.52 ...
## $ TA_max : num [1:5114] -2.14 -3.26 -2.93 6.94 9.91 ...
## $ RH_avg : num [1:5114] 88.5 86.5 87.9 72.6 63.1 ...
## $ SW_IN_avg : num [1:5114] 125 128 145 150 150 ...
## $ NETRAD_avg : num [1:5114] -34.28 -33.93 -8.28 -7.74 3.35 ...
## $ PPFD_IN_avg : num [1:5114] 243 249 282 291 291 ...
## $ PPFD_IN_sum : num [1:5114] 11407 11948 13517 13951 13985 ...
## $ VPD_avg : num [1:5114] 0.0375 0.0569 0.0477 0.1822 0.3185 ...
## $ VPD_max : num [1:5114] 0.226 0.267 0.174 0.543 0.774 ...
## $ LE_avg : num [1:5114] 13.28 8.77 9.01 15.69 28.8 ...
## $ H_avg : num [1:5114] 15.94 7 10.08 2.15 -6.45 ...nrow(seg)## [1] 5114A simple example to start - Average Daily Temperature:
seg %>%
ggplot(aes(x = date, y = TA_avg)) +
geom_line(size = 0.4, color = "tan") +
labs(title = "US-Seg Flux Tower - Average Daily Temperature")
summary(seg)## date GPP_int GPP_dayint RECO_int
## Min. :2007-01-01 Min. :0.0000 Min. :0.00000 Min. :0.001662
## 1st Qu.:2010-07-02 1st Qu.:0.1340 1st Qu.:0.09379 1st Qu.:0.195740
## Median :2013-12-31 Median :0.2662 Median :0.22001 Median :0.345905
## Mean :2013-12-31 Mean :0.5864 Mean :0.53488 Mean :0.521193
## 3rd Qu.:2017-07-01 3rd Qu.:0.7426 3rd Qu.:0.68267 3rd Qu.:0.658740
## Max. :2020-12-31 Max. :5.2217 Max. :5.10939 Max. :5.882377
##
## RECO_dayint RECO_nightint NEE_int NEE_dayint
## Min. :0.0000 Min. :0.00000 Min. :-3.008381 Min. :-3.78218
## 1st Qu.:0.1021 1st Qu.:0.07787 1st Qu.:-0.195228 1st Qu.:-0.34236
## Median :0.1836 Median :0.15992 Median :-0.006409 Median :-0.07767
## Mean :0.2997 Mean :0.22151 Mean :-0.065162 Mean :-0.23520
## 3rd Qu.:0.3883 3rd Qu.:0.27416 3rd Qu.: 0.141722 3rd Qu.: 0.04146
## Max. :4.3728 Max. :1.60844 Max. : 5.810892 Max. : 4.31067
##
## NEE_nightint ET_int ET_dayint P_int
## Min. :-0.85730 Min. :0.01599 Min. :0.0000 Min. : 0.0000
## 1st Qu.: 0.03406 1st Qu.:0.19703 1st Qu.:0.1950 1st Qu.: 0.0000
## Median : 0.11403 Median :0.38021 Median :0.3622 Median : 0.0000
## Mean : 0.17004 Mean :0.58996 Mean :0.5605 Mean : 0.6319
## 3rd Qu.: 0.22711 3rd Qu.:0.75909 3rd Qu.:0.7085 3rd Qu.: 0.0000
## Max. : 1.56843 Max. :4.26213 Max. :4.2112 Max. :46.3441
##
## TA_avg TA_max RH_avg SW_IN_avg
## Min. :-18.914 Min. :-16.44 Min. : 6.062 Min. : 9.236
## 1st Qu.: 6.075 1st Qu.: 13.47 1st Qu.:26.113 1st Qu.:163.145
## Median : 14.267 Median : 21.77 Median :38.656 Median :235.172
## Mean : 13.712 Mean : 20.88 Mean :40.155 Mean :231.808
## 3rd Qu.: 21.993 3rd Qu.: 28.99 3rd Qu.:52.572 3rd Qu.:300.307
## Max. : 30.313 Max. : 38.58 Max. :98.833 Max. :385.824
##
## NETRAD_avg PPFD_IN_avg PPFD_IN_sum VPD_avg
## Min. :-86.76 Min. : 0.0535 Min. : 0 Min. :0.006617
## 1st Qu.: 40.06 1st Qu.: 303.3559 1st Qu.:14138 1st Qu.:0.554238
## Median : 77.19 Median : 442.0772 Median :20836 Median :1.106678
## Mean : 77.84 Mean : 438.4796 Mean :20492 Mean :1.272664
## 3rd Qu.:112.63 3rd Qu.: 576.5788 3rd Qu.:27422 3rd Qu.:1.830674
## Max. :654.11 Max. :1977.1261 Max. :34865 Max. :4.143947
## NA's :195 NA's :100
## VPD_max LE_avg H_avg
## Min. :0.01723 Min. : 0.5317 Min. :-30.02
## 1st Qu.:1.18398 1st Qu.: 6.3996 1st Qu.: 30.39
## Median :2.16125 Median : 12.0348 Median : 50.88
## Mean :2.33556 Mean : 18.5688 Mean : 52.94
## 3rd Qu.:3.34840 3rd Qu.: 23.8515 3rd Qu.: 74.22
## Max. :6.43246 Max. :134.6514 Max. :143.23
## Note that there is no missing data in this data set.
Time period of the US-Seg data:
seg %>%
summarize(Start_Date = min(date),
End_Date = max(date))## # A tibble: 1 × 2
## Start_Date End_Date
## <date> <date>
## 1 2007-01-01 2020-12-31Compute the wavelet power spectrum for Average Temperature time series:
seg_ta_avg.w <- analyze.wavelet(seg, "TA_avg",
loess.span = 0.75,
dt = 1, dj = 1/20,
make.pval = TRUE, method = "white.noise",
n.sim = 25) # should set n.sim higher - default is 100## Smoothing the time series...
## Starting wavelet transformation...
## ... and simulations...
##
|
| | 0%
|
|=== | 4%
|
|====== | 8%
|
|======== | 12%
|
|=========== | 16%
|
|============== | 20%
|
|================= | 24%
|
|==================== | 28%
|
|====================== | 32%
|
|========================= | 36%
|
|============================ | 40%
|
|=============================== | 44%
|
|================================== | 48%
|
|==================================== | 52%
|
|======================================= | 56%
|
|========================================== | 60%
|
|============================================= | 64%
|
|================================================ | 68%
|
|================================================== | 72%
|
|===================================================== | 76%
|
|======================================================== | 80%
|
|=========================================================== | 84%
|
|============================================================== | 88%
|
|================================================================ | 92%
|
|=================================================================== | 96%
|
|======================================================================| 100%
## Class attributes are accessible through following names:
## series loess.span dt dj Wave Phase Ampl Power Power.avg Power.pval Power.avg.pval Ridge Period Scale nc nr coi.1 coi.2 axis.1 axis.2 date.format date.tz# See the WaveletComp analyze.wavelet documentation for more argument options, such as lowerPeriod, upperPeriod, method, etc. Where do you think the most powerful periodicity will be for average temperature?
# you can define your continuous color palette
# mypalette <- "terrain.colors(100)"
wt.image(seg_ta_avg.w,
col.contour = "white",
# color.palette = mypalette,
lwd = 1,
show.date = TRUE,
main = "US-Seg Flux Tower Daily Mean Air Temperature (C)",
legend.params = list(width = 1.8, label.digits = 4))
Average Power of the Series:
# look at the average power of the series
maximum.level = 1.001*max(seg_ta_avg.w$Power.avg) # See WaveletComp Guided Tour for explanation
wt.avg(seg_ta_avg.w, maximum.level = maximum.level,
legend.coords = "bottomright",
main = "Average Power for US-Seg Mean Daily Air Temperature")
Reconstruct the original time series for periods whose average power was found to be significant:
reconstruct(seg_ta_avg.w,
show.date = TRUE,
legend.coords = "bottomright",
ylim = c(-40, 20),
main = "Reconstruction of US-Seg Flux Tower \nDaily Mean Airt Temperature",
col = c('tan', 'black'))## Your input object class is 'analyze.wavelet'...
## Your time series 'TA_avg' will be reconstructed...
## Starting the reconstruction process...
## Original (detrended) and reconstructed series are being plotted...
## Class attributes are accessible through following names:
## series rec.waves loess.span lvl only.coi only.sig siglvl only.ridge rnum.used rescale dt dj Period Scale nc nr axis.1 axis.2 date.format date.tzThe Period information can be used to find and examine periodicities of interest:
seg_ta_avg.w$Period## [1] 2.000000 2.070530 2.143547 2.219139 2.297397 2.378414
## [7] 2.462289 2.549121 2.639016 2.732081 2.828427 2.928171
## [13] 3.031433 3.138336 3.249010 3.363586 3.482202 3.605002
## [19] 3.732132 3.863745 4.000000 4.141060 4.287094 4.438278
## [25] 4.594793 4.756828 4.924578 5.098243 5.278032 5.464161
## [31] 5.656854 5.856343 6.062866 6.276673 6.498019 6.727171
## [37] 6.964405 7.210004 7.464264 7.727491 8.000000 8.282119
## [43] 8.574188 8.876556 9.189587 9.513657 9.849155 10.196485
## [49] 10.556063 10.928322 11.313708 11.712686 12.125733 12.553346
## [55] 12.996038 13.454343 13.928809 14.420007 14.928528 15.454981
## [61] 16.000000 16.564239 17.148375 17.753112 18.379174 19.027314
## [67] 19.698311 20.392970 21.112127 21.856644 22.627417 23.425371
## [73] 24.251465 25.106691 25.992077 26.908685 27.857618 28.840015
## [79] 29.857056 30.909963 32.000000 33.128478 34.296751 35.506223
## [85] 36.758347 38.054628 39.396621 40.785940 42.224253 43.713288
## [91] 45.254834 46.850742 48.502930 50.213382 51.984153 53.817371
## [97] 55.715236 57.680030 59.714111 61.819925 64.000000 66.256955
## [103] 68.593502 71.012446 73.516695 76.109255 78.793242 81.571880
## [109] 84.448506 87.426576 90.509668 93.701485 97.005860 100.426765
## [115] 103.968307 107.634741 111.430472 115.360059 119.428223 123.639850
## [121] 128.000000 132.513910 137.187003 142.024892 147.033389 152.218511
## [127] 157.586485 163.143760 168.897013 174.853153 181.019336 187.402969
## [133] 194.011721 200.853529 207.936613 215.269482 222.860944 230.720118
## [139] 238.856446 247.279700 256.000000 265.027821 274.374006 284.049785
## [145] 294.066779 304.437021 315.172970 326.287521 337.794025 349.706306
## [151] 362.038672 374.805938 388.023441 401.707058 415.873227 430.538965
## [157] 445.721888 461.440237 477.712892 494.559400 512.000000 530.055641
## [163] 548.748013 568.099570 588.133558 608.874043 630.345940 652.575041
## [169] 675.588050 699.412611 724.077344 749.611876 776.046882 803.414116
## [175] 831.746454 861.077929 891.443777 922.880474 955.425783 989.118801
## [181] 1024.000000 1060.111282 1097.496026 1136.199139 1176.267116 1217.748086
## [187] 1260.691879 1305.150082 1351.176101 1398.825223 1448.154688 1499.223753
## [193] 1552.093764 1606.828232 1663.492908The 365-day periodicity is between records 151 and 152. This looks at the power values for periods between 256 and 512 days - octaves 8 and 9.
First, convert power data to a matrix
seg_ta_avg.w_power <- as.matrix(seg_ta_avg.w$Power)The dimensions of the Power matrix have the number of periods as columns and the number of records in the time series as rows
dim(seg_ta_avg.w$Power)## [1] 195 5114It is then possible to pull out various results from analyze.wavelet. For example, looking at the periodicities from 256 to 512:
seg_256_to_512_power <- as.matrix(seg_ta_avg.w_power[141:161,])
dimnames(seg_256_to_512_power) <- list(paste0("period_", round(seg_ta_avg.w$Period[141:161], 2)))
seg_256_to_512_power_t <- as_tibble(t(seg_256_to_512_power))
seg_256_to_512_power_df <- tibble(date = seg$date, seg_256_to_512_power_t)
seg_256_to_512_power_df_long <- seg_256_to_512_power_df %>%
pivot_longer(contains("period_"))
seg_256_to_512_power_df_long %>%
ggplot(aes(x = date, y = value, color = name)) +
geom_line() +
facet_wrap(~ name) +
theme(legend.position="none") +
labs(title = "Average daily temperature periods around 365",
x = "Date",
y = "Power")
Significant p-values:
seg_ta_avg.w_psig <- as.matrix(seg_ta_avg.w$Power.pval)
seg_256_to_512_psig <- as.matrix(seg_ta_avg.w_psig[141:161,])
dimnames(seg_256_to_512_psig) <- list(paste0("period_", round(seg_ta_avg.w$Period[141:161], 2)))
seg_256_to_512_psig_t <- as_tibble(t(seg_256_to_512_psig))
seg_256_to_512_psig_df <- tibble(date = seg$date, seg_256_to_512_psig_t)
seg_256_to_512_psig_df_long <- seg_256_to_512_psig_df %>%
pivot_longer(contains("period_"))
seg_256_to_512_psig_df_long %>%
ggplot(aes(x = date, y = value, color = name)) +
geom_line() +
facet_wrap(~ name) +
theme(legend.position="none") +
labs(title = "Average daily temperature power p-values around 365\n(y-axis free)",
x = "Date",
y = "p-value") +
facet_wrap(~name, scales = "free_y") 
The high power event in 2011
seg %>%
filter(year(date) == 2011) %>%
ggplot(aes(x = date, y = TA_avg)) +
geom_line(size = 0.6) +
labs(title = "2011 Daily Average Temperature (C)")
(start_2011 <- ymd("2011-01-01") - ymd(min(seg$date)))## Time difference of 1461 days(end_2011 <- ymd("2011-12-31") - ymd(min(seg$date)))## Time difference of 1825 daysLook at power bands for periods 16-32
seg_2011_power <- as.matrix(seg_ta_avg.w_power[61:81, 1461:1825])
dimnames(seg_2011_power) <- list(paste0("period_", round(seg_ta_avg.w$Period[61:81], 2)))
seg_2011_power_t <- as_tibble(t(seg_2011_power))
seg_2011_power_df <- tibble(date = seq(from = ymd("2011-01-01"),
to = ymd("2011-12-31"),
by = "1 day"),
seg_2011_power_t)
seg_2011_power_df_long <- seg_2011_power_df %>%
pivot_longer(contains("period_"))
seg_2011_power_df_long %>%
ggplot(aes(x = date, y = value, color = name)) +
geom_line() +
facet_wrap(~ name) +
theme(legend.position="none") +
labs(title = "Average daily temperature in 2011",
x = "Date",
y = "Power")
That was to just show the basic analytic process and how to identify areas of interest in your wavelet images and data.
Another synthetic example that shows how wavelets can detect changes in periodicity.
This example is directly from the guided tour.
Cross-wavelet - analyze two time series variables
Daily Precipitation (mm) - P_int Net Ecosystem Exchange (gC/m^2) - NEE_int
seg %>%
ggplot(aes(x = date, y = P_int)) +
geom_line(size = 0.4, color = "tan") +
labs(title = "US-Seg - Daily Precipitation",
x = "Date",
y = "Precipitation (mm)")
seg %>%
ggplot(aes(x = date, y = NEE_int)) +
geom_line(size = 0.4, color = "tan") +
labs(title = "US-Seg - Net Ecosystem Exchange",
x = "Date",
y = "NEE (gC/m^2)")
seg_p_int.w <- analyze.wavelet(seg, "P_int",
loess.span = 0.75,
dt = 1, dj = 1/20,
make.pval = TRUE, method = "white.noise",
n.sim = 25)## Smoothing the time series...
## Starting wavelet transformation...
## ... and simulations...
##
|
| | 0%
|
|=== | 4%
|
|====== | 8%
|
|======== | 12%
|
|=========== | 16%
|
|============== | 20%
|
|================= | 24%
|
|==================== | 28%
|
|====================== | 32%
|
|========================= | 36%
|
|============================ | 40%
|
|=============================== | 44%
|
|================================== | 48%
|
|==================================== | 52%
|
|======================================= | 56%
|
|========================================== | 60%
|
|============================================= | 64%
|
|================================================ | 68%
|
|================================================== | 72%
|
|===================================================== | 76%
|
|======================================================== | 80%
|
|=========================================================== | 84%
|
|============================================================== | 88%
|
|================================================================ | 92%
|
|=================================================================== | 96%
|
|======================================================================| 100%
## Class attributes are accessible through following names:
## series loess.span dt dj Wave Phase Ampl Power Power.avg Power.pval Power.avg.pval Ridge Period Scale nc nr coi.1 coi.2 axis.1 axis.2 date.format date.tzmypalette <- "terrain.colors(100)"
wt.image(seg_p_int.w,
col.contour = "white",
color.palette = mypalette,
lwd = 1,
show.date = TRUE,
main = "US-Seg Flux Tower Daily Precipitation (mm)",
legend.params = list(width = 1.8, label.digits = 4))
seg_nee_int.w <- analyze.wavelet(seg, "NEE_int",
loess.span = 0.75,
dt = 1, dj = 1/20,
make.pval = TRUE, method = "white.noise",
n.sim = 25)## Smoothing the time series...
## Starting wavelet transformation...
## ... and simulations...
##
|
| | 0%
|
|=== | 4%
|
|====== | 8%
|
|======== | 12%
|
|=========== | 16%
|
|============== | 20%
|
|================= | 24%
|
|==================== | 28%
|
|====================== | 32%
|
|========================= | 36%
|
|============================ | 40%
|
|=============================== | 44%
|
|================================== | 48%
|
|==================================== | 52%
|
|======================================= | 56%
|
|========================================== | 60%
|
|============================================= | 64%
|
|================================================ | 68%
|
|================================================== | 72%
|
|===================================================== | 76%
|
|======================================================== | 80%
|
|=========================================================== | 84%
|
|============================================================== | 88%
|
|================================================================ | 92%
|
|=================================================================== | 96%
|
|======================================================================| 100%
## Class attributes are accessible through following names:
## series loess.span dt dj Wave Phase Ampl Power Power.avg Power.pval Power.avg.pval Ridge Period Scale nc nr coi.1 coi.2 axis.1 axis.2 date.format date.tzwt.image(seg_nee_int.w,
col.contour = "white",
color.palette = mypalette,
lwd = 1,
show.date = TRUE,
main = "US-Seg Flux Tower Net Ecosystem Exchange (gC/m^2)",
legend.params = list(width = 1.8, label.digits = 4))
Cross-wavelet analysis of precipitation and NEE
seg_precip_nee <- seg %>%
select(date, P_int, NEE_int) %>%
as.data.frame()
seg_precip_nee.wc <- analyze.coherency(seg_precip_nee, my.pair = c("P_int", "NEE_int"),
loess.span = 0.75,
dt = 1, dj = 1/20,
make.pval = TRUE,
n.sim = 20)## Smoothing the time series...
## Starting wavelet transformation and coherency computation...
## ... and simulations...## Warning in sqrt(sPower.x * sPower.y): NaNs produced##
|
| | 0%
|
|==== | 5%
|
|======= | 10%
|
|========== | 15%
|
|============== | 20%
|
|================== | 25%
|
|===================== | 30%
|
|======================== | 35%
|
|============================ | 40%
|
|================================ | 45%
|
|=================================== | 50%
|
|====================================== | 55%
|
|========================================== | 60%
|
|============================================== | 65%
|
|================================================= | 70%
|
|==================================================== | 75%
|
|======================================================== | 80%
|
|============================================================ | 85%
|
|=============================================================== | 90%
|
|================================================================== | 95%
|
|======================================================================| 100%
## Class attributes are accessible through following names:
## series loess.span dt dj Wave.xy Angle sWave.xy sAngle Power.xy Power.xy.avg Power.xy.pval Power.xy.avg.pval Coherency Coherence Coherence.avg Coherence.pval Coherence.avg.pval Wave.x Wave.y Phase.x Phase.y Ampl.x Ampl.y Power.x Power.y Power.x.avg Power.y.avg Power.x.pval Power.y.pval Power.x.avg.pval Power.y.avg.pval sPower.x sPower.y Ridge.xy Ridge.co Ridge.x Ridge.y Period Scale nc nr coi.1 coi.2 axis.1 axis.2 date.format date.tzwc.image(seg_precip_nee.wc,
main = "Cross-wavelet power spectrum for US-Seg Precipitation and NEE",
legend.params = list(lab = "cross-wavelet power levels"),
periodlab = "period (days)",
col.contour = "white",
color.palette = mypalette,
lwd = .8,
show.date = TRUE,
col.arrow = "red",
which.image = "wp"
)
wc.sel.phases(seg_precip_nee.wc, sel.period = 192,
only.sig = TRUE,
which.sig = "wc",
siglvl = 0.05,
phaselim = c(-pi,+pi), ## default if legend.horiz = FALSE
legend.coords = "topright", legend.horiz = FALSE,
main = "", sub = "", timelab = "")